## replication code based on canada_v4.r, timestamp 20201019.
## this code is very slow and "base R" and hard to extend (not how I would do now). but it works.
# Fundamentally the issue is that we need to generate tau and the outcome measure separately for each imputation and each way in which we can measure preferences, beliefs, and vote choice.  
# In the code below we do this on the fly, cycling over preference measures (mclogit-based, just party, pre vs post), beliefs (forecast vs results) and vote choice measures (pre vs post); in the same loop we run regressions and combine results across imputations. 

# I considerably sped up the analysis by only doing a baseline situation plus one-step deviations from it rather than all permutations. 

# If I were to do again, it would be a lot tidier and it would have a lot of purrr::map stuff applied to imputations. 

### see data_prep_imputation_part.r for dataset construction:
  ### preferences drawn from optimal weighting of party, leader, candidate, w imputation
  ### preferences drawn from party only, w imputation
  ### preferences drawn from party only, no imputation 
  ### preferences drawn from post-election party ratings
  ### preferences drawn from post-election leader ratings
### this script does the analysis for each dataset (set of imputations, really), combines the results, outputs figures. 

output.date = "20201019" # 20201027" # 20191027" # 20190321"
library(lfe)
source("util_vote_value_functions.r")
source("util_plotting_functions_v2.r")
source("util_combine_imputation_models.r")

load(file = "lists_for_analysis.RData")

# exclude the small number of ridings where "other" got significant votes 
results = read.csv("vote_shares_by_riding_2015.csv")
ridings_to_exclude = results$constituency_number[results$oth > .05]  # this is 5 constituencies, 473 respondents. 


parties = c("lib", "con", "ndp", "bq", "grn")
ss = c(10, 25, 40, 53, 70, 85)
pchs = c(19, 21, 17)
this.ymax = .7

pp.mat.list = list()

for(s in ss){
  cat("s = ", s, ": storing pivotal probabilities.\n", sep = "")
  pp.mat.list[[paste0("s=", s)]] = list("results" = list(), "forecast" = list())
  pp.mat.results = read.csv(paste0("pivot_probs_x1M_2015_canada_s_", s, "_results.csv"))
  pp.mat.forecast = read.csv(paste0("pivot_probs_x1M_2015_canada_s_", s, ".csv"))
  for(i in 1:nrow(pp.mat.results)){
    pp.mat.list[[paste0("s=", s)]][["results"]][[paste0("c", pp.mat.results[i,1])]] = convert_pp_vec_to_pp_mat(pp.mat.results[i,-1])
    pp.mat.list[[paste0("s=", s)]][["forecast"]][[paste0("c", pp.mat.results[i,1])]] = convert_pp_vec_to_pp_mat(pp.mat.forecast[i,-1])
  }
}


## idea is to get a single big for loop to do the analysis, rather than a series in which we store and move on. 
## maybe extract some functions. 

# big.r.list = responsiveness.reg.list = list()
binary.offsets = c(-.05, .05)
ternary.offsets = c(-.075, 0, .075)

belief.sources = c("results", "forecast") 
preference.types = names(lists.for.analysis) # includes pes party and pes leader # also "weighted.utilities"  "party.utilities"     "raw.party.utilities"

noc = 10 # number of categories of tau

intx.sto = responsiveness.sto = list()
sample.size.list = list()
reg.list.sto = list()
tau.vals.sto = list()

sr.mat.sto = list()
data.store = list()

# these are the baseline values: we do these and one-step deviations from these
baseline_s <- 53
baseline_outcome.suffix <- "both"
baseline_preference.type <- "weighted.utilities"
baseline_belief.source <- "results"

for(belief.source in belief.sources){
  for(preference.type in preference.types){
    for(s in ss){
      
      # for efficiency, restrict to cases we show in the paper: we show baseline and one-step variations from it. 
      variations_from_baseline <- c(s != baseline_s, preference.type != baseline_preference.type, belief.source != baseline_belief.source)
      if(sum(variations_from_baseline) > 1){next}
      
      
      cat(" === ", belief.source, "-based beliefs; ", preference.type, "-based preferences; s = ", s, " === \n", sep = "")
      this.pp.mat.list = pp.mat.list[[paste0("s=", s)]][[belief.source]]
      # we need to store the results so we can combine across imputations 
      leaf.list = list("age" = list(), "edu" = list(), "income" = list(), "gender" = list(), "leftright" = list())
      intx.reg.list = responsiveness.reg.list = intx.reg.list.no.tau.cats = intx.reg.list.attention = intx.reg.list.cpr = list("both" = leaf.list, "pre" = leaf.list, "post" = leaf.list)

      for(imputation.number in 1:length(lists.for.analysis[[preference.type]])){
        
        for(outcome.suffix in c("both", "pre", "post")){responsiveness.reg.list[[outcome.suffix]][[imputation.number]] = list()}
        
        # get this imputation for this dataset.
        d = lists.for.analysis[[preference.type]][[imputation.number]]
        
        if(imputation.number == 1 & belief.source == belief.sources[1] & preference.type == preference.types[1] & s == ss[1]){
          write(paste0("Before analysis we have ", nrow(d), " observations.\n"), file = "missingness_log.txt", append = T)
        }
        
        # can't vote BQ outside Quebec
        d$bq_therm[d$region != "Quebec"] = NA
        d$vote_party_post[d$region != "Quebec" & d$vote_party_post == "bq"] = NA
        d$vote_party_both[d$region != "Quebec" & d$vote_party_both == "bq"] = NA
        
        # this actually doesn't change across imputations
        party_id_mat = vote_mat_both = vote_mat_cps = vote_mat_pes = NULL
        for(party in parties){
          party_id_mat = cbind(party_id_mat, as.integer(d$party_id == party))
          vote_mat_both = cbind(vote_mat_both, as.integer(d$vote_party_both == party))
          vote_mat_cps = cbind(vote_mat_cps, as.integer(d$vote_party_pre == party))
          vote_mat_pes = cbind(vote_mat_pes, as.integer(d$vote_party_post == party))
        }
        
        # this does vary across imputations, because the therm scores are imputed. 
        u_mat_raw = d[,paste0(parties, "_therm")]
        u_mat = u_mat_raw + party_id_mat/100  # just a bit of a push for the party_id party to break any ties 
        colnames(u_mat) = parties
        
        fave_party_mat = matrix(NA, nrow = nrow(u_mat), ncol = length(parties))
        row.max = apply(u_mat, 1, max, na.rm = T)
        for(j in 1:ncol(u_mat)){
          fave_party_mat[,j] = as.integer(u_mat[,j] == row.max)
        }
        
        # so does this because it depends on therm scores
        d$leftright = NA
        d$leftright[u_mat[,which(parties == "lib")] > u_mat[,which(parties == "con")]] = 1
        d$leftright[u_mat[,which(parties == "lib")] < u_mat[,which(parties == "con")]] = 0
        
        d$left.fp = NA
        d$left.fp[fave_party_mat[, which(parties == "lib")] == 1] = 2
        d$left.fp[fave_party_mat[, which(parties == "ndp")] == 1] = 1
        # d$left.fp[fave_party_mat[, which(parties == "grn")] == 1] = 3
        
        # some exclusions 
        vote_recorded_both = apply(vote_mat_both[,1:5], 1, sum, na.rm = T) == 1  
        number_of_faves = apply(fave_party_mat, 1, sum, na.rm = T)
        max_u_vec = apply(u_mat, 1, max, na.rm = T)
        
        ### identify rows to include 
        use = !is.na(max_u_vec) & !d$constituencynumber %in% ridings_to_exclude & number_of_faves == 1 & vote_recorded_both ### continue here if necessary
        
        if(imputation.number == 1 & belief.source == belief.sources[1] & preference.type == preference.types[1] & s == ss[1]){
          write(paste0(sum(is.na(max_u_vec)), " have no max_u_vec -- probably no non-missing ratings.\n"), file = paste0("missingness_log", output.date, ".txt"), append = T)
          write(paste0(sum(d$constituencynumber %in% ridings_to_exclude), " are in a riding with more than 5% support for a minor party.\n"), file = paste0("missingness_log", output.date, ".txt"), append = T)
          write(paste0(sum(number_of_faves == 0), " don't have a favorite party.\n"), file = paste0("missingness_log", output.date, ".txt"), append = T)
          write(paste0(sum(number_of_faves > 1), " have more than one favorite party.\n"), file = paste0("missingness_log", output.date, ".txt"), append = T)
          write(paste0(sum(!vote_recorded_both), " do not have a vote recorded in either PES or CPS.\n"), file = paste0("missingness_log", output.date, ".txt"), append = T)
          write(paste0("Combining these, ", sum(!use), " are excluded at this stage and ", sum(use), " are included.\n"), file = paste0("missingness_log", output.date, ".txt"), append = T)
        }
        
        # define argument list for responsiveness regression: we cycle through this to make our figures
        argument.list = list(
          "all" = list(cat.mat = matrix(use, ncol = 1), cat.labels = "All", offsets = c(0)),
          "gender" = list(cat.mat = cbind(use & d$gender == 2, use & d$gender == 1), cat.labels = c("Women", "Men"), offsets = binary.offsets),
          "income" = list(cat.mat = cbind(use & d$income_cat == 3, use & d$income_cat == 2, use & d$income_cat == 1), cat.labels = c("High income", "Medium income", "Low income"), offsets = ternary.offsets),
          "edu" = list(cat.mat = cbind(use & d$edu_cat == 3, use & d$edu_cat == 2, use & d$edu_cat == 1), cat.labels = c("University plus", "Some higher education", "High school or less"), offsets = ternary.offsets),
          "age" = list(cat.mat = cbind(use & d$age_cat == 3, use & d$age_cat == 2, use & d$age_cat == 1), cat.labels = c("Age 56+", "Age 37-56", "Age 36 or below"), offsets = ternary.offsets),
          "left.fp" = list(cat.mat = cbind(use & d$left.fp == 1, use & d$left.fp == 2), cat.labels = c("First pref: NDP", "First pref: Liberal"), offsets = ternary.offsets)
        )
        
        #### generate tau, vote codes, etc. 
        vv.mat = matrix(NA, ncol = length(parties), nrow = nrow(d))
        tau.vec = rep(NA, nrow(d))
        vote.code.both = vote.code.pre = vote.code.post = rep(NA, nrow(d))
        
        k = 0  # counter
        cat("Imputation ", imputation.number, " -- filling in ", sum(use), " rows: ", sep = "")
        for(i in which(use)){
          k = k + 1
          if(k%%1000 == 0){cat(k, " ")}
          out = get_vote_values_and_tau_and_vote_code_for_voter_given_utilities_and_pivotal_prob_mat_and_vote_indicator_vec(u_mat[i,], this.pp.mat.list[[paste0("c", d$constituencynumber[i])]], vote_mat_both[i,1:5])
          vv.mat[i,] = out$vv
          tau.vec[i] = out$tau
          vote.code.both[i] = out$vote_code
          vote.code.post[i] = vote_code_from_utilities_and_vote_indicator_vec_and_vote_value(u_mat[i,], vote_mat_pes[i,1:5], out$vv)
          vote.code.pre[i] = vote_code_from_utilities_and_vote_indicator_vec_and_vote_value(u_mat[i,], vote_mat_cps[i,1:5], out$vv)
        }
        
        # tau cats: depend on tau
        decile.cuts = quantile(tau.vec[!is.na(tau.vec) & use], probs = seq(0, 1,.1))
        decile.cuts[abs(decile.cuts) == min(abs(decile.cuts))] = 0
        tau.cats = NA
        decile.midpoints = c()
        for(i in 1:(length(decile.cuts) - 1)){
          tau.cats[use & tau.vec >= decile.cuts[i] & tau.vec < decile.cuts[i+1]] = i
          decile.midpoints = c(decile.midpoints, mean(decile.cuts[i:(i+1)]))
        }
        
        # assemble dataset for analysis
        E = data.frame(tau.vec = tau.vec, tau.cats = as.factor(tau.cats), constituencynumber = d$constituencynumber, age_cat = as.factor(d$age_cat), income_cat = as.factor(d$income_cat), edu_cat = as.factor(d$edu_cat), gender = as.factor(d$gender), left.fp = as.factor(d$left.fp), attention = d$attention_1, cpr = d$correctly_predicted_result)
        
        ### Store it here to avoid having to rerun imputations every time?
        data.store[[paste0(belief.source, "_", preference.type, "_s_", s, "_imp_", imputation.number)]] = E
        
        cat("Done.\nRunning regressions and storing regressions for this imputation: ")
      
        for(outcome.suffix in c("both", "pre", "post")){
          
          # define the outcome for this run 
          E$tactical.vote = I(get(paste0("vote.code.", outcome.suffix)) == 2)
          
          this_name <- paste0(belief.source, "_", preference.type, "_s_", s, "_", outcome.suffix)
          if(!this_name %in% names(sr.mat.sto)){
            sr.mat.sto[[this_name]] = list()
          }

          ## put that matrix in this list
          this_mat <- matrix(NA, nrow = 2, ncol = 2)
          this_mat[1,] <- c(sum(use & E$tau.vec <= 0), sum(use & E$tau.vec > 0))
          this_mat[2,] <- c(sum(use & E$tau.vec <= 0 & E$tactical.vote), sum(use & E$tau.vec > 0 & E$tactical.vote))
          sr.mat.sto[[this_name]][[imputation.number]] = this_mat

          if(imputation.number == 1 & s == ss[1]){
            sample.size.list[[paste0(preference.type, "_", outcome.suffix, "_age_cat")]] = sum(use & !is.na(E$tau.vec) & !is.na(E$tactical.vote) & !is.na(E$age_cat))
            sample.size.list[[paste0(preference.type, "_", outcome.suffix, "_edu_cat")]] = sum(use & !is.na(E$tau.vec) & !is.na(E$tactical.vote) & !is.na(E$edu_cat))
            sample.size.list[[paste0(preference.type, "_", outcome.suffix, "_income_cat")]] = sum(use & !is.na(E$tau.vec) & !is.na(E$tactical.vote) & !is.na(E$income_cat))
            sample.size.list[[paste0(preference.type, "_", outcome.suffix, "_gender")]] = sum(use & !is.na(E$tau.vec) & !is.na(E$tactical.vote) & !is.na(E$gender))
            sample.size.list[[paste0(preference.type, "_", outcome.suffix, "_left_fp")]] = sum(use & !is.na(E$tau.vec) & !is.na(E$tactical.vote) & !is.na(E$left.fp))
          }

          # run interaction regressions, which I will combine via Rubin rules below to make a single figure
          intx.reg.list[[outcome.suffix]][["all"]][[imputation.number]] = felm(tactical.vote ~ tau.cats + I(tau.vec > 0) | 0 | 0 | constituencynumber, data = E[use,])
          intx.reg.list[[outcome.suffix]][["age"]][[imputation.number]] = felm(tactical.vote ~ tau.cats + I(tau.vec > 0)*age_cat | 0 | 0 | constituencynumber, data = E[use,])
          intx.reg.list[[outcome.suffix]][["income"]][[imputation.number]] = felm(tactical.vote ~ tau.cats + I(tau.vec > 0)*income_cat | 0 | 0 | constituencynumber, data = E[use,])
          intx.reg.list[[outcome.suffix]][["edu"]][[imputation.number]] = felm(tactical.vote ~ tau.cats + I(tau.vec > 0)*edu_cat | 0 | 0 | constituencynumber, data = E[use,])
          intx.reg.list[[outcome.suffix]][["gender"]][[imputation.number]] = felm(tactical.vote ~ tau.cats + I(tau.vec > 0)*gender | 0 | 0 | constituencynumber, data = E[use,])
          intx.reg.list[[outcome.suffix]][["left.fp"]][[imputation.number]] = felm(tactical.vote ~ tau.cats + I(tau.vec > 0)*left.fp | 0 | 0 | constituencynumber, data = E[use,])
          
          # same thing but with no flexible control for tau
          intx.reg.list.no.tau.cats[[outcome.suffix]][["all"]][[imputation.number]] = felm(tactical.vote ~ I(tau.vec > 0) | 0 | 0 | constituencynumber, data = E[use,])
          intx.reg.list.no.tau.cats[[outcome.suffix]][["age"]][[imputation.number]] = felm(tactical.vote ~ I(tau.vec > 0)*age_cat | 0 | 0 | constituencynumber, data = E[use,])
          intx.reg.list.no.tau.cats[[outcome.suffix]][["income"]][[imputation.number]] = felm(tactical.vote ~ I(tau.vec > 0)*income_cat | 0 | 0 | constituencynumber, data = E[use,])
          intx.reg.list.no.tau.cats[[outcome.suffix]][["edu"]][[imputation.number]] = felm(tactical.vote ~ I(tau.vec > 0)*edu_cat | 0 | 0 | constituencynumber, data = E[use,])
          intx.reg.list.no.tau.cats[[outcome.suffix]][["gender"]][[imputation.number]] = felm(tactical.vote ~ I(tau.vec > 0)*gender | 0 | 0 | constituencynumber, data = E[use,])
          intx.reg.list.no.tau.cats[[outcome.suffix]][["left.fp"]][[imputation.number]] = felm(tactical.vote ~ I(tau.vec > 0)*left.fp | 0 | 0 | constituencynumber, data = E[use,])
          
          # controlling for attention
          intx.reg.list.attention[[outcome.suffix]][["all"]][[imputation.number]] = felm(tactical.vote ~ tau.cats + I(tau.vec > 0)*attention | 0 | 0 | constituencynumber, data = E[use,])
          intx.reg.list.attention[[outcome.suffix]][["age"]][[imputation.number]] = felm(tactical.vote ~ tau.cats + I(tau.vec > 0)*age_cat + I(tau.vec > 0)*I(attention > 7) | 0 | 0 | constituencynumber, data = E[use,])
          intx.reg.list.attention[[outcome.suffix]][["income"]][[imputation.number]] = felm(tactical.vote ~ tau.cats + I(tau.vec > 0)*income_cat + I(tau.vec > 0)*I(attention > 7) | 0 | 0 | constituencynumber, data = E[use,])
          intx.reg.list.attention[[outcome.suffix]][["edu"]][[imputation.number]] = felm(tactical.vote ~ tau.cats + I(tau.vec > 0)*edu_cat + I(tau.vec > 0)*I(attention > 7) | 0 | 0 | constituencynumber, data = E[use,])
          intx.reg.list.attention[[outcome.suffix]][["gender"]][[imputation.number]] = felm(tactical.vote ~ tau.cats + I(tau.vec > 0)*gender + I(tau.vec > 0)*I(attention > 7) | 0 | 0 | constituencynumber, data = E[use,])
          intx.reg.list.attention[[outcome.suffix]][["left.fp"]][[imputation.number]] = felm(tactical.vote ~ tau.cats + I(tau.vec > 0)*left.fp + I(tau.vec > 0)*I(attention > 7) | 0 | 0 | constituencynumber, data = E[use,])
          
          # controlling for correctly predicted result
          intx.reg.list.cpr[[outcome.suffix]][["all"]][[imputation.number]] = felm(tactical.vote ~ tau.cats + I(tau.vec > 0)*cpr | 0 | 0 | constituencynumber, data = E[use,])
          intx.reg.list.cpr[[outcome.suffix]][["age"]][[imputation.number]] = felm(tactical.vote ~ tau.cats + I(tau.vec > 0)*age_cat + I(tau.vec > 0)*cpr | 0 | 0 | constituencynumber, data = E[use,])
          intx.reg.list.cpr[[outcome.suffix]][["income"]][[imputation.number]] = felm(tactical.vote ~ tau.cats + I(tau.vec > 0)*income_cat + I(tau.vec > 0)*cpr | 0 | 0 | constituencynumber, data = E[use,])
          intx.reg.list.cpr[[outcome.suffix]][["edu"]][[imputation.number]] = felm(tactical.vote ~ tau.cats + I(tau.vec > 0)*edu_cat + I(tau.vec > 0)*cpr | 0 | 0 | constituencynumber, data = E[use,])
          intx.reg.list.cpr[[outcome.suffix]][["gender"]][[imputation.number]] = felm(tactical.vote ~ tau.cats + I(tau.vec > 0)*gender + I(tau.vec > 0)*cpr | 0 | 0 | constituencynumber, data = E[use,])
          intx.reg.list.cpr[[outcome.suffix]][["left.fp"]][[imputation.number]] = felm(tactical.vote ~ tau.cats + I(tau.vec > 0)*left.fp + I(tau.vec > 0)*cpr | 0 | 0 | constituencynumber, data = E[use,])
          
          # now the responsiveness.regs, which are used to make the SRF figures
          # for each imputation number and each outcome suffix and each element of argument list above I want a list of regressions
          for(element in names(argument.list)){
            use.mat = argument.list[[element]]$cat.mat
            this.list = list()
            responsiveness.reg.list[[outcome.suffix]][[imputation.number]][[element]] = list()
            for(y in 1:ncol(use.mat)){
              responsiveness.reg.list[[outcome.suffix]][[imputation.number]][[element]][[y]] = felm(tactical.vote ~ tau.cats - 1 | 0 | 0 | constituencynumber, data = E[use.mat[,y],])
            }
          }
        }  # cycling across outcomes
        cat("Done.\n")
      }  # cycling across imputations 
      cat("done.\nCombining regression output and storing/printing figures:\n")
      # for intx we have to store it because the figure combines results for different preference types.  
      for(outcome.suffix in c("both", "pre", "post")){
        # intx first
        this.intx.reg.list = this.intx.reg.list.no.tau.cats = this.intx.reg.list.attention = this.intx.reg.list.cpr = this.reg.list.sto = leaf.list
        for(characteristic in c("all", "age", "edu", "income", "gender", "left.fp")){
          this.intx.reg.list[[characteristic]] = combine_imputation_models(intx.reg.list[[outcome.suffix]][[characteristic]])
          this.intx.reg.list.no.tau.cats[[characteristic]] = combine_imputation_models(intx.reg.list.no.tau.cats[[outcome.suffix]][[characteristic]])
          this.intx.reg.list.attention[[characteristic]] = combine_imputation_models(intx.reg.list.attention[[outcome.suffix]][[characteristic]])
          this.intx.reg.list.cpr[[characteristic]] = combine_imputation_models(intx.reg.list.cpr[[outcome.suffix]][[characteristic]])
        }
        # for intx, we store 
        intx.sto[[paste0(belief.source, "_", preference.type, "_", s, "_", outcome.suffix)]] = this.intx.reg.list 
        intx.sto[[paste0(belief.source, "_", preference.type, "_", s, "_", outcome.suffix, "_no_tau_cats")]] = this.intx.reg.list.no.tau.cats 
        intx.sto[[paste0(belief.source, "_", preference.type, "_", s, "_", outcome.suffix, "_attention")]] = this.intx.reg.list.attention 
        intx.sto[[paste0(belief.source, "_", preference.type, "_", s, "_", outcome.suffix, "_cpr")]] = this.intx.reg.list.cpr 
        
        
        # for responsiveness we can make the plot here. 
        # the tricky thing is that we need to combine the elements across imputations, but imputations are not the lowest level of the tree.
        # really I should output the imputations and operate from there. 
        # but it would still be tricky if I did that. 
        # for each characteristic, we need a list with 2 or 3 elements (depending on how many categories), each of which is matrix with term, estimate, std.error columns. and we plot that.
        for(characteristic in c("all", "age", "edu", "income", "gender", "left.fp")){
          reg.list.to.plot = list()
          for(category_index in 1:length(responsiveness.reg.list[[outcome.suffix]][[1]][[characteristic]])){
            regs.to.combine = list()
            for(imputation.index in 1:length(responsiveness.reg.list[[outcome.suffix]])){
              regs.to.combine[[imputation.index]] = responsiveness.reg.list[[outcome.suffix]][[imputation.index]][[characteristic]][[category_index]]
            }
            responsiveness.sto[[paste0(characteristic, "_", s, "_", preference.type, "_", outcome.suffix)]] = regs.to.combine
            reg.list.to.plot[[category_index]] = combine_imputation_models(regs.to.combine)
          }
          
          # for replication, this is one long filename
          pdf(paste0(belief.source, "_responsiveness_", "group_comparison_by_", gsub("\\.", "_", characteristic), "_s_", s, "_saturated_10_bins_", gsub("\\.", "_", preference.type), "_", outcome.suffix, ".pdf"), width = 3.5, height = 3)
          par(mfrow = c(1,1), mar = c(4,5,1,1), cex = .9)
          
          make.plot.of.saturated.regressions.2(reg.list.to.plot, plot.confints = T, error.bar.confint.plot = F, x.lab = "", x.labels = rep("", noc), pchs = pchs, tones = c(.4, .6, .8), offsets = c(0, 0, 0), y.lims = c(0, this.ymax), y.lab = "Proportion voting best insincere")
          mtext("Tactical incentive (10 bins, even spacing)", 1, line = 1.75, cex = .8)
          if(length(argument.list[[characteristic]]$cat.labels) > 1){legend("topleft", legend = argument.list[[characteristic]]$cat.labels, pch = pchs[1:length(argument.list[[characteristic]]$cat.labels)], col = "black", text.col = "black", lty = 1:length(reg.list.to.plot), cex = .8)}
          
          dev.off()
          
          this.reg.list.sto[[characteristic]] = reg.list.to.plot
          
        }
        key <- paste0(belief.source, "_s_", s, "_", gsub("\\.", "_", preference.type), "_", outcome.suffix)
        tau.vals.sto[[key]] = tapply(tau.vec, tau.cats, mean)
        reg.list.sto[[key]] = this.reg.list.sto
        
      }
      cat("Done.\n\n")
    }
  }
}

# for all ss and beliefs and outcomes, comparing across 3 models (weighted, party ratings, party ratings no imp) 
for(s in baseline_s){
  for(belief.source in baseline_belief.source){
    for(outcome.suffix in c("both", "pre", "post")){
      reg.list = list()
      for(preference.type in preference.types[1:3]){
        reg.list[[preference.type]] = intx.sto[[paste0(belief.source, "_", preference.type, "_", s, "_", outcome.suffix)]]
      }
      pdf(paste0(belief.source, "_by_preference_type_all_characteristics_s_", s, "_", outcome.suffix, ".pdf"), width = 9, height = 5.5)
      par(mar = c(4,5,1,1))
      make_intx_plot(reg.list)
      dev.off()      
    }
  }
}

# for all ss and beliefs and outcomes, comparing across all preference types -- see how this looks, consider for main paper.   
for(s in baseline_s){
  for(belief.source in baseline_belief.source){
    for(outcome.suffix in c("both", "pre", "post")){
      reg.list = list()
      for(preference.type in preference.types){
        reg.list[[preference.type]] = intx.sto[[paste0(belief.source, "_", preference.type, "_", s, "_", outcome.suffix)]]
      }
      pdf(paste0(belief.source, "_by_preference_type_including_pre_post_all_characteristics_s_", s, "_", outcome.suffix, ".pdf"), width = 9, height = 5.5)
      par(mar = c(4,5,1,1))
      make_intx_plot(reg.list, model.labels = c("Model 1: utility from party, leader, and candidate ratings", "Model 2: utility from party ratings only", "Model 3: utility from party ratings only (no imputation)", "Model 4: utility from post-election party ratings", "Model 5: utility from post-election leader ratings"), offsets = seq(-.2,.2, length = 5), pchs = c(19,21,2:4), ylims = c(-.1, .3), legend.loc = "topleft")
      dev.off()      
    }
  }
}


# just pre and post prefs 
for(s in baseline_s){
  for(belief.source in baseline_belief.source){
    for(outcome.suffix in c("both", "pre", "post")){
      reg.list = list()
      for(preference.type in preference.types[c(1,4,5)]){
        reg.list[[preference.type]] = intx.sto[[paste0(belief.source, "_", preference.type, "_", s, "_", outcome.suffix)]]
      }
      pdf(paste0(belief.source, "_by_preference_type_only_pre_post_all_characteristics_s_", s, "_", outcome.suffix, ".pdf"), width = 9, height = 5.5)
      par(mar = c(4,5,1,1))
      make_intx_plot(reg.list, model.labels = c("Model 1: utility from pre-election party, leader, and candidate ratings", "Model 2: utility from post-election party ratings", "Model 3: utility from post-election leader ratings"), offsets = seq(-.2,.2, length = 3), pchs = c(19,21,2), ylims = c(-.1, .35), legend.loc = "topleft")
      dev.off()      
    }
  }
}

## assemble and output strategic responsiveness matrix 
mat <- sr.mat.sto[["results_weighted.utilities_s_53_both"]][[1]]
for(j in 2:5){
  mat <- mat + sr.mat.sto[["results_weighted.utilities_s_53_both"]][[j]]
}
mat <- mat/5
mat <- rbind(mat, mat[2,]/mat[1,])
mat[1:2,] <- round(mat[1:2,], 0)
mat[3,] <- round(mat[3,], 3)
rownames(mat) <- c("Number of observations", "Number casting best insincere vote", "Proportion casting best insincere vote")
s <- rep(NA, nrow(mat))
for(j in 1:nrow(mat)){
  s[j] <- paste(c(rownames(mat)[j], mat[j,]), collapse = "&")
}
write(paste(s, collapse = " \\\\ \n"), file = "sr_matrix.tex")


# for one s and model and beliefs, comparing across outcomes 
for(s in baseline_s){
  for(preference.type in baseline_preference.type){
    for(belief.source in baseline_belief.source){
      reg.list = list()
      for(outcome.suffix in c("both", "pre", "post")){
        reg.list[[outcome.suffix]] = intx.sto[[paste0(belief.source, "_", preference.type, "_", s, "_", outcome.suffix)]]
      }
      pdf(paste0(belief.source, "_by_outcome_suffix_all_characteristics_s_", s, "_", gsub("\\.", "_", preference.type), ".pdf"), width = 9, height = 5.5)
      par(mar = c(4,5,1,1))
      make_intx_plot(reg.list, model.labels = c("Vote choice from post-election or pre-election survey", "Vote choice from pre-election survey", "Vote choice from post-election survey"), legend.loc = "topright", ylims = c(-.1, .2))
      dev.off()      
    }
  }
}

# across ss 
for(outcome.suffix in c("both", "pre", "post")){
  for(preference.type in baseline_preference.type){
    for(belief.source in baseline_belief.source){
      reg.list = list()
      for(s in ss){
        reg.list[[paste0("s=", s)]] = intx.sto[[paste0(belief.source, "_", preference.type, "_", s, "_", outcome.suffix)]]
      }
      pdf(paste0(belief.source, "_by_s_all_characteristics_", outcome.suffix, ".pdf"), width = 9, height = 5.5)
      par(mar = c(4,5,1,1))
      make_intx_plot(reg.list, model.labels = paste0("s=", ss), offsets = seq(-.3, .3, length = length(ss)), pchs = c(6, 2:3, 19, 4:5), horiz = T)
      dev.off()      
    }
  }
}

# across belief types 
for(s in baseline_s){
  for(preference.type in baseline_preference.type){
    for(outcome.suffix in c("both")){
      reg.list = list()
      for(belief.source in c("results", "forecast")){
        reg.list[[belief.source]] = intx.sto[[paste0(belief.source, "_", preference.type, "_", s, "_", outcome.suffix)]]
      }
      pdf(paste0("by_belief_source_all_characteristics_s_", s, "_", gsub("\\.", "_", preference.type), "_", outcome.suffix, ".pdf"), width = 9, height = 5.5)
      par(mar = c(4,5,1,1))
      make_intx_plot(reg.list, model.labels = c("Beliefs centered on results", "Beliefs centered on forecasts"), offsets = c(-.1, .1), ylims = c(-.15, .2))
      dev.off()      
      cat("??")
    }
  }
}

# a plot showing results when different controls are added 

controls = c("attention", "cpr", "no_tau_cats")
for(outcome.suffix in "both"){
  for(preference.type in baseline_preference.type){
    for(belief.source in baseline_belief.source){
      reg.list = list()
      reg.list[["none"]] = intx.sto[[paste0(belief.source, "_", preference.type, "_", s, "_", outcome.suffix)]]
      for(control in controls){
        reg.list[[control]] = intx.sto[[paste0(belief.source, "_", preference.type, "_", s, "_", outcome.suffix, "_", control)]]
      }
      pdf(paste0(belief.source, "_by_control_variable_all_characteristics_", outcome.suffix, ".pdf"), width = 9, height = 5.5)
      par(mar = c(4,5,1,1))
      make_intx_plot(reg.list, model.labels = c("Baseline",  "With control for pol. attention", "With control for pol. knowledge", "Without quasi-deciles of tau"), offsets = seq(-.3, .3, length = length(controls) + 1), pchs = c(19, 1:(length(controls))), horiz = F)
      dev.off()      
    }
  }
}



